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We have formulated a kinetic theory for a condensed atomic gas in a trap, i.e., a generalized 
Gross-Pitaevskii equation, as well as a quantum-Boltzmann equation for the normal and anomalous 
fluctuations [R. Walser et al., Phys. Rev. A, 59, 3878 (1999)]. In this article, the theory is applied 
to the case of an isotropic configuration and we present numerical and analytical results for the 
reversible real-time propagation, as well as irreversible evolution towards equilibrium. 
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I. INTRODUCTION 

More than 70 years ago, S. Bose and A. Einstein pro- 
posed a provocative hypothesis — that at ultra-low tem- 
peratures a novel state of matter should exist. They pre- 
dicted this state could be attained by cooling an ordinary 
gas towards absolute zero. At a well-defined point in 
this process, a spontaneous transition should occur and 
change the state of matter from an unordered ensemble of 
individual particles into one collective entity. This single 
object, now devoid of its many-particle character, ought 
to evolve as a collective matter wave. 

With the discovery of superfluidity in liquid helium in 
1938 and its subsequent explanation in terms of Bose- 
Einstein condensation (BEC), the hypothesis had been 
firmly established. In turn, this phenomenon has had a 
major impact on the development of modern quantum 
physics. Today, BEC is fundamental to our understand- 
ing of many low-temperature phenomena and it is the 
cornerstone of many quantitative explanations. However, 
up to 1995, condensation of a weakly interacting, atomic 
Bose-Einstein gas had never been achieved, as such. 

With the ground-breaking accomplishment of condens- 
ing atomic ^^Rb by E. Cornell and C. Wieman et al. [|l[ , of 
sodium by W. Ketterle et al. ||] , and lithium by R. Hulet 
et al. 1^] , a new chapter of quantum statistical physics has 
been opened. For the first time, it is possible to study 
in a table-top experiment quantum statistical effects of 
material objects on a human scale (up to 5 mm) — the 
very phenomena that govern the otherwise microscopic 
physics of nuclear matter, macroscopic quantum liquids, 
or astronomical objects, such as neutron stars. 

Today, more bosonic alkali elements have crossed the 
transition temperature, in particular atomic hydrogen 
as well as ®^Rb and many more vastly improved ex- 
periments have been carried out. For example, it is now 
possible to examine multi-component condensates [^0, 
to create vortices and to prepare novel topological 

modes For a list of current experiments see Ref. 

PI , or the review article in Ref. . However, the tech- 
nological breakthrough of combining laser cooling with 



evaporative cooling is not limited to bosonic species only. 
Most recently, the fermionic isotope of potassium ""^K has 
also been cooled successfully below the Fermi tempera- 
ture [|l|. 

Instigated by these spectacular experiments, strongly 
renewed interest has developed in their quantitative de- 
scription. While cold quantum gases had been studied 
extensively in the 1950-60s, they were mainly considered 
as precursor theories for strongly interacting systems, 
such as liquid helium. Thus, most of the available results 
were focused on spatially uniform systems in thermal 
equilibrium. Excellent accounts of these standard results 
can be found, for example, in the textbooks and mono- 
graphs [p"^-p9|. However, the spatial non-uniformity, the 
thermal isolation resulting from the confinement in a 
ultra-high vacuum trap, as well as the large disparity 
of collision and relaxation time-scales, are indispensable 
ingredients for a quantitative description of today's ex- 
periments. 

To account for these differences that distinguish the 
present experimental situation from the homogeneous 
Bose-Einstein gas [20 -23 , a growing number of equilib- 
rium and nonequilibrium kinetic theories have been re- 
cently presented ||2^-[29|]. However, the effort to go be- 
yond the mean-field description of the Gross-Pitaevskii 
equation is considerable. Thus, the research for a 
unified description of the equilibrium and nonequilibrium 
situation is still very active. 

In this article, we explore numerically and analytically 
some of the implications of the reversible and irreversible 
evolution of a condensed gas immersed in the noncon- 
densate cloud. The points discussed are organized as 
follows. Section || revisits the main results of our kinetic 
theory |3l[ |, i. e., the two-particle Hamiltonian and the en- 
ergy and number conserving coUisional kinetic equations 
for the condensate, as well as the normal and anoma- 
lous fluctuations. In Sec. Ill, we specialize these kinetic 



equations for a completely isotropic situation. Based on 
these prerequisites, we discuss in Sec. IV the results of 



propagating the coUisionless mean-field and the Hartree- 
Fock-Bogoliubov (HFB) equations in real-time. Finally, 
in Sec. we study the evolution of an ergodic distribu- 
tion towards equilibrium in the presence of collisions. 
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II. KINETIC MASTER EQUATIONS 



B. Dynamical evolution 



A. Master variables 

The kinetic master equation of the weakly interact- 
ing dilute atomic gas describes the coupled evolution of 
the condensed fraction immersed in the quantum fluctu- 
ations. In this context, we associate the condensate with 
a c- number field a^(t) that represents the expectation 
value of the quantum field {a^(t)). The field operator 
removes a particle from point x and satisfies the scalar, 
equal-time commutation relation, 

[a^,al] ^ 5{x-y), (f) 

of a boson. The position representation {|x)} used above 
is not necessarily the most suitable basis to formulate a 
kinetic theory. It proves to be more useful to postpone 
the choice of a particular representation and to formulate 
the theory in terms of a general single particle basis {|ji)} 
that spans the same single-particle Hilbert space: 

«x = E%H*i)- (2) 

il 

In the case of an unstructured (scalar) atomic conden- 
sate, three external quantum labels (zi) are sufficient to 
describe its motional state in space, completely.^ In this 
manner, we can expand any field as 

(a) = a^^ In) EE ai |1) ee a. (3) 

il 

Here we have simplified the notation by dropping the 
name of the dummy variable, i. e., ii = 1, and by assum- 
ing implicit summation over repeated indices, as usual. 

In an analogous fashion, we can describe the normal 
density of the atomic gas / — (a^a) = /'^'^^ -I- / by a 
Hcrmitian tensor operator of rank (1,1): 

/ = /i2 |1)(2|, /W^a^ai |1)(2|. (4) 

Moreover, we will always decompose any quantum av- 
erage into a mean-field contribution and the remaining 
fluctuations. Similarly, we define the anomalous aver- 
ages m = (aa) = m*-*^-* -|- m as symmetric tensors of rank 
(2,0), 

m = mi2 |1)|2), m(^) = a2ai |1) |2) , (5) 
and their symmetric conjugates as n = to*2 (1| (2|. 



*This is readily generalized to accommodate multiple inter- 
nal electronic configurations if ii encompasses more quantum 
labels accordingly, i. e., — \ni, h, mi; Fi, Ah, . . .). 



The kinetic evolution of a weakly interacting gas is 
primarily governed by the motion of the individual par- 
ticles in the external trapping potential and by binary 
collisions. Simultaneous collisions of more than two par- 
ticles are unlikely events in a dilute gas. Consequently, 
we will disregard such processes and use the following 
number-conserving Hamiltonian operator: 

H = ij (") + if = a\a^ + 01234 a\ala^a^. (6) 

Here, -ff'"-* denotes a single-particle Hamiltonian op- 
erator with matrix elements iJ^*^-* — (I|p2/(2r7i) + 
Vcxt (x) 1 2) . For the external trapping potential, we as- 
sume a three-dimensional isotropic harmonic oscillator, 
T4xt(x) = mw^ (x^ + + 2:^)72. In most of the present 
experiments with large, stable condensates, the two-body 
interaction potentials Vbin(xi — X2) are repulsive and of 
short range. From such potentials, we can obtain two- 
particle matrix elements as 

01234 ^ l(^) ^ Fbin(xi - X2) |3) ® |4) , (7) 
^1234 ^^1243 ^^2134^^2143^ (8) 

Only the symmetric part of the two-particle matrix el- 
ement 01234 jg physically relevant. Therefore, we have 
explicitly (S) symmetrized it. In the low kinetic energy 
range that we are interested in, s-wave scattering is the 
dominant two-particle scattering event ]3^|3^. Thus, by 
discarding all details of the two-particle potential, we can 
describe the interaction strength with a single parameter 
Vq related to the scattering length aghy Vq — Airh'^as/m. 
This limit corresponds to a singular interaction potential, 
i. e., Vbin(xi, X2) — Vb ^(xi — X2). In the case of this delta 
potential, one finds for the two-body matrix elements: 

^1234 ^ j rf3^(i|^)(2|x)(x|3)(x|4), (9) 

which need not be symmetrized, as they are symmetric 
already. However, considering the caveats that are re- 
lated to the singular functional form of the two-particle 
potential p5[ | , we will only rely on the existence and sym- 
metry of the two-particle matrix elements as defined in 
Eq. (0). 

C. Mean field equations 

Based on these assumptions, we have derived a set of 
kinetic equations that describe the dynamical evolution 
of the condensate fraction immersed in a cloud of non- 
condensate particles. By discarding all of the interactions 
except for the condensate's self-interaction, they reduce 
to the familiar Gross-Pitaevskii (GP) equation for the 
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mean field a. However, due to the presence of anoma- 
lous fluctuations to, this nonlinear, but otherwise unitary 
GP equation acquires a contribution proportional to the 
time-reversed or complex conjugated field a*. 

To represent these equations compactly, it is useful to 
arrange them in a 2 x 2 matrix form. Moreover, we trans- 
form this field equation to a frame co-rotating with a 
positive frequency fi defined by a{t) = exp (— i/it) 
However, in order not to overload the notation, we will 
suppress the overlinc in the following generalized GP- 
equation 



dt 



X 



(10) 



The two-component state vector x — (o^iO;*)^? intro- 
duced above keeps track of the forward and time-reversed 
components of the mean-field. It is symmetric under 
time-reversal, i.e., x — '^iX* ■ The Pauli matrix tri 
achieves the exchange of upper and lower components 
and is defined in Appendix ^ . 

Two distinct processes govern the real-time evolution 
of the mean-field. First, there is the generalized GP- 
propagator that is defined as 



H = 



(11) 



The two contributions that define this symplectic propa- 
gator are a normal Hermitian Hamiltonian operator 

n^ = i/("'+lC/^,e)-H2[/;-/X, (12) 

as well as a symmetric anomalous coupling strength 

Ua = V^. (13) 

It is easy to identify Hf^ with the well known unitary 
GP-propagator that accounts for the free evolution of the 
mean- field (i?*^"^ — /i), its self-interaction C//(c) , as well as 
the energy shift Uj caused by the presence of the non- 
condensate cloud. However, due to the existence of the 
anomalous fluctuations there is also a coupling through 
to the time-reversed field. For convenience, we have 
introduced two auxiliary operators Uf and Vm ■ Explic- 
itly, they are defined in terms of the two-body matrix 
elements, such as 



C/; = 2 0i2'3'4'/3'2' |1) (4'|, 
and a first order anomalous coupling strength 

V,^^2^^^'^'^'m,,,, |1)|2'). 



(14) 



(15) 



Second, there are all of the coUisional second-order damp- 
ing rates and energy shifts |20-2^] that are given by 



T< 



(16) 



and the time-reversed contribution = — criT<*iTi. It 
can be shown that they are equivalent to the extended 
Beliaev rates |Q. 

The forward and backward transition rates T^, T^, 
T^, and T^, describe the bosonically enhanced scatter- 
ing of noncondensate particles into and out of the con- 
densate. In turn, these transition rates are formed from 
various binary scattering processes F, and are given by 



^A^ - ff{l+f) 



2 ^ffhn ' 



^X- - r(i+/)(i+/)/ ■ 



2r 



(l + /)mn' 



and 



T< 
T> 



2r/s(i+/), 



2F 



(17) 

(18) 

(19) 
(20) 



Within the Born-Markov approximation of kinetic the- 
ory, we define these elemental collision processes as 





^12'3'4' 




'2' 


'3' 


''"/3'l"/4'2"/4"2' |1>(3"|, 


T/m/ = 8 


^12'3'4' 




'2' 


'3' 


'*73'l"m4.3"/4"2' |1) |2") 




^12'3'4' 




'2' 


'3' 


'*"/3'l"m4'3"n2"2' |1) (4"| 


^mmn — ^ 


^12'3'4' 




'2' 


'3' 


TO3/4//TO4'3//7l2"2' 1) l" 



(21) 

During a binary collision event, two particles can con- 
serve their energy only approximately. After all, the in- 
dividual scattering event happens within a medium and 
the asymptotics can not be reached within the finite du- 
ration of the collision. Thus, within the limits of the 
Born-Markov approximation, any second order collision 
operator accrues a dispersive as well as a dissipative part 
from the complex valued matrix-element 



/1"2"3"4" 



i,l"2"3"4" 



1 



7]- Z Al"2' 



3"4" 



(22) 



It is essentially non-zero only if the energy difference 
Ai"2"3"4" = El" {t) + £2" (t) — £3" {t) — £4" (t) between the 
pre- and post-collision energies is smaller than an energy 
uncertainty 77: 



lim — 

»)^o+ 7] — i A 



^5^(A) + zn,^ 



(23) 



On general physical grounds, it can be argued that this 
uncertainty rj is bracketed by the binary collision rate, on 
one side, and the energy uncertainty arising from the fi- 
nite duration of an individual collision event on the other 
side. As we have shown, one has also the liberty to choose 
a more accurate intermediate propagator such that the 
single particle energies e {t) and the eigen-states incorpo- 
rate mean field shifts. 
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D. Normal and anomalous fluctuations 

The normal and anomalous fluctuations f{t) and m (t) 
of a quantum field are not independent quantities, but 
actually they are the components a generalized single- 
time density operator G-* (t) : 



G> = 



I 



[i + fY 



> 0. 



(24) 



The non-negativity of this co-variance operator implies 
that the magnitude of the anomalous fluctuations is lim- 
ited by the normal depletion through a Cauchy-Schwartz 
inequality (see Appendix ^). In the general context 
of Green function's Jl5| , p^ , this single-time density op- 
erator (t) can also be viewed as a particular limit 
of a time-ordered (T), two-time Green function G{T,t), 
i.e., G^{t) = liuiT^t^T G{T,t). Consequently, it is 
also necessary to consider the opposite limit and to de- 
fine a time-reversed, single-time density operator through 
G^(t) = limr-tt_ T G{T,t). Explicitly, this operator is 
given by 



ai G> Vi = G= 



as 



l + f m 
n f* 



(25) 



With the help of these definitions, we can now present the 
results of the kinetic theory as a generalized Boltzmann 
equation for the single-time density operator G^ {t) as 



dt 



G' 



-iY.G' 



r< G< - r> G- 



h. c. 



(26) 



In analogy to the previous discussion of the mean-field 
dynamics, we again find that the evolution of the density 
operator is ruled by two types of propagators. First, there 
is the Hartree-Fock-Bogoliubov (HFB) self-energy oper- 
ator S that can be obtained also by variational methods 
p6[ . In detail, this symplectic self-energy is given by 



E = 



Ea/- E^ 



(27) 



where we have introduced Hermitian Hamiltonian oper- 
ators and symmetric anomalous coupling potentials as 



— ^(m(<=)+rfi)- 



(28) 
(29) 



It is important to note the different weighing factors of 
the mean- field potential in Eqs. ( |l2|) and ([28|), as well 
as the appearance of the anomalous condensate density 
m'-'^^ in Eq. (p9|). This HFB-operator is the usual starting 
point of any finite-temperature calculations. Depending 
on additional considerations, i. e ., 'ga. pless vs. conserv- 
ing approximations' (see Refs. [^2|j3^^-|40[] ) , the anoma- 
lous couplings Vffi are usually discarded from Eqs. ( p^ ) 
and (29). However, since we do go beyond a first or- 
der calculation, we need to retain all contributions for 
consistency. 



Second, the Boltzmann equation, Eq. (|26|), introduces 
forward and backward collision operators F*" and F^. 
They are responsible for particle transfer out of and into 
the condensate on one hand, and lead to thermal equi- 
libration within the non-condensate cloud, on the other 
hand. These forward and backward collision operator are 
defined by 



F< = 



^ A/- 

-F>* 



r< 
r> * 



(30) 



and F^ = — criF^*cri, where 
F 



- )/(!+/) 



2 F 



(/+/(<^))mn 



F 



(l+/+/(c))(i+/)/ 1 (i+/)/(c)/ 1 



F 



(31) 



2 F 



{l+f+f<-'=y)fhn 



(1 + /)?TI<=)?1 



(l+/)mri(' 



, (32) 



and 



__ip_ __l'P__ 



^ '^(i+/+/<'^')m/ + ^(i+/)™(<=)/ + r(i+/)r7i/("^) ) ■ (34) 



It is interesting to note that all of the collision processes 
that contribute to the Boltzmann equation, Eq. (|2^) , are 
of the same basic structure as the collision operators in 
the GP equation, Eq. ([lo|). In particular, one can gener- 
ate all of the processes F< and F-* by functional differen- 
tiation from T< and T> . This very fact is actually is the 
key principle to the functional-analytic Green functions 
method described in Ref. ||l5| and, for example, leads to 
the gapless Beliaev approximation | |2^ , p| . 



E. Conservation laws 

1. Number 

The total particle number A'^ is a conserved quantity if 
the atoms evolve under the generic two-particle Hamilto- 
nian operator H given by Eq. (^, i. e., [H, N] = 0. This 
conservation law implies that the system is invariant un- 
der a global phase change a ^ a exp {i ip) . By using 
this continuous symmetry, i. e., a — > a exp (iip), f f, 
and m ^ fh exp (2 i (p), it is easy to see that the kinetic 
equations Eqs. (^0|) and ( p^ ) are also explicitly number 
conserving at all times: 



{N{t)) = Tr{/(^)(t)} + Tr{/(i)} = const. 



(35) 



Nevertheless, it is important to note that there are al- 
ways coherent and incoherent processes present that do 
transfer particles between the condensate and the non- 
condensate clouds, continuously. 
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2. Energy 

In the absence of any time-dependent external driv- 
ing fields, such as optical lasers or magnetic rf- fields, 
the overall energy _ff must be conserved as well. To 
find the expectation value of the total system energy 
E = = Tr{HfT(i)}, one can use the same power- 
series expansion of the coarse-grained many-particle den- 
sity matrix a{€) that leads to the kinetic equations. Thus, 
within the limits of the Born-Markov approximation and 
the systematic application of Wick's theorem, we have 
obtained first and second order contributions for the en- 
ergy E = Tr{H {af^^^^^^ +af^^^^^)} + 0[3]. Explicitly, this 

energy functional E = E'-"^ + E[f\ + E[ m], is given as 

= Tr{[i?(0) + 1(1 [/^(., + 2 C/; + z (T<. - T>.))] /(^) 

+ \{Vfn+^{T^-T>,))n('''^}, (36) 

E[f] = Tr{[H(") + 1 (2 Uf,., +2U^~ + t (T^ - T^))] /}, 

(37) 
(38) 



E[m] = Tr{- (y„(e,+,^, + * (F^ - T^)) n}. 



For example, the same first order results can be found in 
Ref. [hw, derived by a variational procedure. 



III. A COMPLETELY ISOTROPIC SYSTEM 



the desired radial symmetry. Thus, according to Refs. 
p2| , ^ , we introduce irreducible representations of tensor 
fields of rank: (2,0), (1,1), (0,2) as 



\nilimi) (n2/2"i2| , 



m2) 



(39) 

The quantum labels carry additional overlines or under- 
lines to indicate whether a function depends only on two 
of the three quantum labels, for example, 1 = (ni,Zi) or 
2 = {l2,iTi2). With these definitions, it is easy to verify 
the following orthogonality relationships: 

Tr{T^(i2)r^,(3 4)^} = <5i3 J^i -^n' <5 w , (40) 



Tr 



{s;„(1 2) 5^,(3 4)^} =<5i3 524'5«'<5w (41) 



In the case of a scalar field {I — 0), one can simplify 
Eq. (|3^) by the following relation for the Clebsch-Gordan 
coefficients C^'/^^q = (-l)'i-'"VV2l7TT. 
Provided the basis states transform under a coordinate 
rotation, here denoted by 7?,, according to the finite di- 
mensional representation of the rotation group I?' (7?.) 



Uk \nlm) = J2 1"^"^') ; 



(42) 



In the previous section, we have reviewed the main re- 
sults of the kinetic theory that describes the coupled evo- 
lution of the condensate immersed in the non-condensate. 
The formal derivation did not rely on a particular trap- 
ping geometry, nor a special form for the binary inter- 
action potential. In order to gain deeper understand- 
ing of the intricate interactions, we will now specialize 
the theory to the most simple, though realistic, three- 
dimensional model: a completely isotropic configura- 
tion, a spherically symmetric harmonic trapping poten- 
tial, and a binary short-range s-wave scattering potential. 



it follows that the set of tensors {T^|| \m\ < 1} and 
{5',J| \m\ < 1} are irreducible as well: 

Un ® Un-' 2) = ^ Ti,(f 2) Vl,^{TZ), (43) 

m' 

Un ® Un Sl{l 2) = ^ SUl 2) Vl,^{n). (44) 
B. Isotropic two- particle matrix element 



A. Irreducible tensor fields 

Complete isotropy is easily achieved for the mean-field 
by decomposing it in terms of a few zero angular mo- 
mentum partial waves. For this purpose, we will use a 
set of basis states {|f)} that can be characterized by ra- 
dial and angular momentum quantum numbers (n, /, m), 
i.e., (x|f) =i?r;(r)il\(0,^). 

However, in order to isolate the isotropic components 
of the non-condensate fluctuations, / and m, we need 
to generalize the concept of partial waves and introduce 
an irreducible set of tensor fields. Furthermore, by only 
selecting the scalar component {I — 0), we can enforce 



The most commonly used model for a short-range bi- 
nary interaction potential is the s-wave hard-core delta 
potential Vbin (xi , X2 ) = Vb ^ (xi — X2 ) . This model poten- 
tial is most suited to describe the low energetic collision 
dynamics of two real particles. However, it has to be 
used with caution in connection with infinite summation 
over virtual, high energy states. It is clear that the en- 
ergy independent scattering approximation fails above a 
certain energy range when the spatial scale of variation 
of the high-energetic wave functions begin to sample the 
detailed form of the interaction potential. Thus, the true 
value of the interaction matrix element ought to decrease 
much faster with energy than the value obtained from the 
simple hard core delta-potential approximation. It is well 
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known that indiscriminate use of the energy independent 
approximation leads to a non physical ultra-violet diver- 
gence ||,||||. 

Considering these limitations, we will use the s-wave 
scattering matrix element that is obtained from the en- 
ergy independent approximation, i.e., 



(.1234 



2 



d3x(l|x)(2|x)(x|3)(x|4). 



(45) 



only for energies below a certain level and truncate it 
appropriately otherwise. 

By using the basis states (x|l) = R\'^{r)Yj-,l_^{9,ip), we 
can decompose the matrix element 0^^^"* into a reduced 
radial part 0^^^^ and a purely geometric factor Y----, 



1,1234 



1234 yl 2 3 4 



(46) 



where 



f,i234 ^ Y2. 
2 



r'drRl^{r)Rl^{r)R^l^ir)Rl-{r) (47) 



and 



00 /-^/l/2^/^^3^4^/^^1^2^ ^1^1 

000 000 m\m2 (mi +m2 ) 77137714 (7711+7712) 

/-O 



47r(2; + l)[nti(2^. + l)]-^/^ 



(48) 



In the context of evaluating the collision integrals, it 
will be necessary to consider products of two matrix 
elements summed over all energetically accessible sub- 
levels. Within the isotropic model, all magnetic sub- 
levels are energetically degenerate. Moreover, spherical 
symmetry demands an equal population distribution and 
rules out the existence of coherences within magnetic sub- 
manifolds. Thus, it will be required to know the magni- 
tude of |yi234|2 averaged over all the magnetic quantum 
numbers. For later reference, we will now introduce such 
conveniently scaled factors as 

/2^((2?i + l)(2/2 + l))^/V27r, (49) 
5^2 34 ^2 (^2 ^iV')"' E 



9 ^12 „3 4 l'-"000 '-^OOO ) 



(50) 



These coupling strengths g^^^^ in Eq. (|50| ) mea- 
sure the amount and principle connectivity between 
pre-coUision and post-collision angular momenta sub- 
manifolds (Z3, 14) —^ (/i, 12). In particular, it establishes a 
parity selection rule such that the coefficients are non- 
vanishing only if the sum of the angular momenta is 
^1 + ^2 + ^3 + ^4 = even. In addition, transitions are al- 
lowed only if the angular momentum I is within a range 
of max(|Zi - ^2!, 1^3 - ^4|) < I < min(Zi + ^2, ^3 + ^4)- 



C. Scalar component of states and energies 

With the help of the auxiliary results established in 
the previous section, we are now able perform the de- 
sired multi-pole decomposition of the kinetic equations. 
The postulate of complete isotropy then implies that we 
can focus on the scalar component of the field {I = 0) 
exclusively. This is 

m = mi2^0(i2), n = nr2S°\l2). ^ ' 

Analogously, we can decompose all normal operators, 
such as the bare single particle Hamiltonian operator, 
as 

i/(o) = i/(0)j4,T(?(14'). (52) 

The first order normal mean-field potential U f and the 
anomalous coupling strength Vm are then 

Uf = ^™' g~''5i^,i^,h-,.T°,{ll'), (53) 
= ^3-4' 5o"(l 2'). (54) 

Finally, the normal and anomalous coUisional contribu- 
tions simplify to 

T/// = ^'^'^>r'"'"'" 5^'''''' -^v^.. 

X /3'l"/4'2"/4"2'r°(13"), (55) 
T /I 2'3'4' / l"2"3"4" 12'3'4'r r r 

X /3,i„m4.3"n2.2'r°(14"), (56) 

P _ /1 2'3'4'il"2"3"4'' 12'3'4'r c c 

X /3a"m4,3"/4"2'5°(12"), (57) 

-p il 2'3'4' i,l"2"3"4" 12'3'4'c r r 

^rnmn—V (Pjj 3 61^,1^,, Ol^,l^,, Ol^,,l^, 

X m3,4„m4,3„n2..2'5'S(ll")- (58) 



IV. REVERSIBLE EVOLUTION 

In this section, we will examine several limiting situa- 
tions of the reversible evolution in order to elucidate the 
complex behavior of the condensed gas. Since canonical 
transforms and Hartree-Fock-Bogoliubov (HFB) opera- 
tors are crucial for an understanding of the reversible 
evolution, we will review the main results ||l6|| . Subse- 
quently, we are going to examine the stationary equilib- 
rium, as well as the reversible real-time evolution of the 
condensed gas. 



A. Structure of the generalized density matrix 

The definition of a generalized density matrix G, i. e., 
either or G^, was given in Eqs. (24) and (p5|). Its spe- 
cific structure implies various important physical proper- 
ties. 
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First of all, we have to assume that there is a basis 
that diagonalizes this (2n x 2 n)-dimensional fluctuation 
matrix. Exactly n of its 2 n eigen- values correspond to 
the positive occupation numbers of finding a particle or, 
more generally, a quasi-particle in a certain mode. For 
a given, but otherwise arbitrary, G matrix, one can con- 
struct this basis by studying the transformation law of 
the density matrix under a canonical transformation T 
(see Appendix ^), 

G' = TGT'^. (59) 

It is important to note that this is not the transformation 
law of a general matrix under coordinate change. This 
would require that = . However, by only using 
the properties of the symplectic transformations, one can 
show that a canonical eigen-value problem is defined by 

(^3 G) = ((73 G') . (60) 

The solution of this eigen-value problem yields the eigen- 
vector matrix and the corresponding diagonal eigen- 
value matrix a^G'. Here, we have introduced standard 
Pauli spin matrices cti and as, which exchange upper and 
lower component of a 2 n dimensional vector, or flip the 
sign of the lower segment, respectively. All normalizable 
states can be rescaled such that Ta^T^ = (73. Now, we 
are able to reconstruct the positive G matrix 

G:^VPV\ (61) 

from its eigen- vectors V — aj,T^ and the diagonal, posi- 
tive occupation number matrix P = (73 G'a^. 

Second, an important feature of an admissible fluc- 
tuation matrix is its consistency with the commutation 
relation, i.e., {aid\) = +"^12 and {0,^(12) = (0.2(3.1). 

This can be expressed compactly as 

(7iG*(7i-G = (73. (62) 

By invoking the properties of a unitary symplectic trans- 
formation, one can show that the elements of the diago- 
nal occupation number matrix P are not 2 n independent 
variables. Actually half of them are determined by the 
other half, P(„+i,...,2n) = 1 + ..,«), or 

(7lP(7l -P=(73. (63) 

In other words, by separating the occupation numbers 
P and the eigen-vector matrix V into a first and second 
half, i.e.(P+,l + P+) = P and {V+,V-) = V, one can 
then decompose a general fluctuation matrix as 

G^V^P^Vl + V_{l + P^)Vl. (64) 



B. Structure of the Hartree-Fock-Bogoliubov 
operator 

The symplectic HFB operator arises not only naturally 
in kinetic theories or variational calculations, but in many 



other contexts involving stability analysis. In the case of 
bosonic fields, the self-energy operator is of the generic 
form: 

In here, Y,j\f stands for a Hermitian operator E^v" = ^jv 
and denotes an anomalous coupling term that has to 
be symmetric E^ — Ej^. The relative size of the opera- 
tors Yijsf and E^ determines the character of the energy 
spectrum. It can either be real valued with pairs of pos- 
itive and negative eigen-energies, or one finds a doubly 
degenerate zero eigen-value, if the energy diS'erence be- 
tween the smallest positive and highest negative vanishes 
(gapless spectrum). In the general case, there is a mixed 
spectrum consisting of pairs of real sign-reversed as well 
as pairs of complex conjugated eigen- values. The eigen- 
vectors W are normalizable with respect to the indefinite 
norm \W\^ = W^a^W, except for those that belong to 
zero or complex eigen-values. It is important to note 
that this energy basis W is in general distinct from the 
instantaneous basis V that diagonalizes the fluctuation 
matrix G in Eq. (|6l]). They do coincide only in equilib- 
rium. The mathematical properties of the eigen-states 
W can be derived easily from the intrinsic symmetries of 
the HFB operator: 

E==-criEVi, (66) 
Et=(73E(73. (67) 

Thus, if W and E are the solutions of the right eigen- 
value problem, 

J:W = WE, (68) 

it follows directly from Eq. (^) that W — a\ W* , is 
also a right eigen-vector but corresponds to the eigen- 
value E = —E*. Starting from the second symmetry in 
Eq. (j^ and the right eigen-value problem of Eq. (|6^), 
it is easy to construct the left eigen-vectors W = (73 
that correspond to the eigen-values E = E*: 

WT,^E*W. (69) 

Finally, from a combination of the results for the right 
and left eigen-vectors, it follows that the eigen-vectors 
are orthogonal with respect to the metric a^: 

O^iE* -E')wlasWE,, (70) 

if E* ^ E' . On the other hand, this relation implies also 
that eigen-vectors that belong to complex eigen-values 
must have zero norm. 

The situation of a doubly degenerate zero energy eigen- 
value E — needs special attention. One can view this 
case as a limit when two non-degenerate states approach 
each other. However, as the energy gap decreases, the 
two eigen-states become more and more collinear. Thus, 
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in the limit of a vanishing energy separation, the dimen- 
sion of the spanned vector space collapses from 2 to 1 and 
E becomes defective. In the present context however, we 
did not encounter this situation (see Ref. jl^ for details). 



C. Stationary solution of the 
Hartree-Fock-Bogoliubov equations 

In spite of the complex nonlinear interactions taking 
place within the atomic gas, the kinetic evolution is com- 
pletely reversible if we disregard all collisionally induced 
redistributions of quasi-particles. Thus, for the moment, 
we will eliminate the collision operators T and T from 
the kinetic Eqs. ( [lO| ) and (26) and we will study the col- 
lisionless stationary equilibrium, as well as the real-time 
evolution in this section. 

With these assumptions, we are left with the following 
set of stationary equations for the mean field x Q'Hd tlis 
fluctuations G^: 



G> 



(71) 
(72) 



The self-energies of the condensate 11 and the noncon- 
densate E are nonlinearly coupled and implicitly include 
the rotation frequency of the mean-field /i. 

However, the equilibrium solution to Eqs. ( [7l| ) and ( |7^ ) 
is not fully determined as it stands. From the results 
of the previous section, we know that any fluctuation 
matrix that is diagonal with respect to the positive 
and negative energy eigen- vectors W = (W+, W-.) of E, 
will be a stationary and complete solution of Eq. ([72|) 



G= 



w. p, wi + W_(1 + P,) wl 



(73) 



By choosing a canonical Bose-Einstein distribution 
P{E > Eo) — l/[exp(/3i?) — 1] for the quasi-particles 
above the non-degenerate ground state Eq > Q and a 
vanishing ground state occupation number P{Eq) = 0, 
we obtain a variationally minimal energy solution for the 
total system at some inverse temperature P ]T^ : 

G>^ J2 PiE)WEWl + {l + P{Ej)W_EWlE- (74) 

E>Eo 

In order to understand the self-consistent equilibrium 
solution of Eqs. ( [7l|) and ([72|), it is useful to examine 
first the potentials that govern the evolution of the con- 
densate, as well as the noncondensate. In Fig. |l|, we de- 
pict the potential energy densities of the normal Hamil- 
tonian operators and Ejv" versus radius that arise 
for the zero angular momentum manifold I = 0, i.e., 
Vcxt + !£//(<:) and Kxt + 2Uf(c), respectively. They are 
compared to the bare isotropic harmonic oscillator po- 
tential Voxt = r'^/'^, for reference. Here and in all of 
the subsequent results, we will use the experimental data 



of a typical ®^Rb condensate |^^. All physical param- 
eters are scaled in the natural units for a harmonic os- 
cillator, i.e., the angular frequency uj = 27r200Hz, the 
atomic mass mgr = 86.9092 amu, the ground state size 
a-H — [fi/{i-^ msr)]^^"^ = 763 nm, the s-wave scattering 
length as = 5.82 nm = 7.63 10^^ gh, a very low temper- 
ature of ksT — 0.2 huj, and a condensate number chosen 
as N^'^') ~ 10'*. The isotropic particle densities are nor- 
malized to iV = drr^ f{r). From the effective cou- 
pling parameter k — as /an, one obtains an estimate of 
the mean-field energy shift as ^tf = (15A^'^'^^k)^/^/2 = 
8.36. This gives an excellent approximation of the self- 
consistent chemical potential of /i = 8.52, as can be seen 
in Fig. |l|. 




FIG. 1. Self-consistent potential energy densities of the 
condensate Hamilton operator IXat (solid), the noncondensate 
Hamilton operator Eat (dashed dot), and the bare harmonic os- 
cillator potential (dashes) versus radius. Energy and length are 
scaled in the natural units for a harmonic oscillator. 



The position densities of the condensate f^''\r), as 
well as the local densities of the normal and anoma- 
lous fluctuations, i.e., /(r) and fh{r), are represented 
in Fig. ||. First of all, it has to be noted that the solu- 
tions are real valued. This important fact follows from 
the detailed structure of the stationary Eqs. ( fzi] ) and 
([Z^), which are invariant under a global phase change. 
Second, if we focus on the condensate density, one can 
see that it closely follows the Thomas-Fermi approxi- 
mation /tf(^) — (''tf ~ ''^)/(2'«)i for radii less than 
rxF = \/2/i « 4.12. This limit is valid in the strong 
coupling regime iV^'^' k > 1, where the kinetic energy is 
a negligible contribution compared to the external trap- 
ping potential and the self-energies. The self-consistent 
solutions for the normal and anomalous densities are de- 
picted in the lower half of Fig. |^. While all normal densi- 
ties are necessarily positive, the anomalous fluctuations 
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carry a negative sign. The anomalous fluctuations are 
the response of the noncondensate medium to a phase- 
coherent mean-field. In analogy to the polarization of 
an atom that is subjected to an electric field, it tries to 
compensate for the external perturbation. In the evalu- 
ation of the normal and anomalous fluctuations, we have 
truncated the flnite temperature sums beyond the radial 
and angular momentum quantum numbers = 14 and 
1 — 6. This leads to fully converged values of the normal 
fluctuations. However, it has to be noted that the values 
of the anomalous fluctuations are still subject to change 
(further decrease). As long as one keeps adding vac- 
uum contributions at an unaltered strength of the s-wave 
matrix-element (jf-lhulhulhulu [ggg ^^as would 

lead to the well known ultra-violet divergence ||2^,^,Q . 
Thus, a judicious truncation, or alternatively an energy 
renormalization is needed to remove the non-physical di- 
vergence that arises solely from the energy-independent 
approximation of the scattering amplitudes. 



barrier l{l + \)/r'^ . The corresponding eigen-energies are 
depicted in Fig. ^. In the context of spatially homoge- 
neous condensed matter systems, these eigen- functions 
are associated with particle-like excitations. 
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FIG. 2. Position density of the mean-field f^^'ir) (solid), the 
Thomas-Fermi approximation /|p'(r) (dashed-dot), as well as 
the normal fluctuations f{r) (solid) and the anomalous fluctua- 
tions fh(r) (dashed-dot). Density and length are scaled in the 
natural units for a harmonic oscillator. 



In Fig. ^ we show a few selected radial eigen-functions 
(x|l) = i?" {r)Yj^{d,(p) of the noncondensate Hamilto- 
nian operator (ij(°) + 2Uf(.)) |i) = £i The lowest 
energy state RJ^^^^{r) is localized at the rim of the con- 
densate rxF = M ~ 4.12. It has a smaller spatial 
extent than the condensate and consequently a higher 
energy. All s-wave functions (/ — 0) have a finite value at 
the origin in contrast to the I > states that must be van- 
ishing at r = 0. An eigen-function that is characterized 
by quantum numbers {n^ ,1) has — 1 nodes. Eigen- 
functions corresponding to higher angular momenta I are 
shifted outwards due to the increased angular momentum 



-0.5' ' ' ' ' ' ' ' ' 
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r 

FIG. 3. Radial eigen-functions R" (r) of the particle-like ba- 
sis states |1) associated with the noncondensate Hamiltonian 
operator H'^'^^ + 2U j{a) versus radius. Depicted are a few rep- 
resentative states for the radial and angular quantum numbers 

n'' = 1,2, 10 and l = Q,2. 



The other relevant set of eigen-states arises from the 
condensate Hamiltonian operator, i.e., the stationary 
GP-equation (ifW -M i/^tc)) |l(^)) = ei jl^'^)). The lowest 
self-consistent energy eigen-state defines the condensate 
wave function. A selection of these eigen-states are shown 
in Fig. ^. As these states correspond to the low energetic 
excitation-modes of the condensate they are referred to 
as phonon-like. The eigen-energies of the isotropic (/ — 0) 
modes are shown in Fig. ^ 

We have compiled the four important positive energy 
spectra that arise in the problem in Fig. In essence, 
these are the spectra of the condensate and the noncon- 
densate Hamiltonians, lij^ and Ejv, as well as their gen- 
eralization in terms of the HFB self-energy operators H 
and S. It can be seen that the s-wave energies of Ha/' 
and H are virtually identical. In contrast to this, one 
finds that the excitation frequencies of the fiuctuations 
S are characteristically shifted downwards from the en- 
ergies of the noncondensate Sa^. It is also important to 
note that the self-energy of the fluctuations Ejv' includes 
the energy shifts of the noncondensate itself. These nu- 
merical results compare well within the limits of validity 
with the perturbative and semi-classical approximations 
of Refs. 

The spectrum of eigen-values of S exhibits a charac- 
teristic energy gap above the condensate energy level /i. 
It is well known that this gap energy vanishes asymptot- 
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ically for a homogeneous system in the thermodynamic 
Umit. By deliberately excluding the anomalous coupling 
Vm strength from a first order theory (Popov approxima- 
tion) or by including the second order Beliaev correction 
of Eq. (^6|) in the self-energies, one can obtain a gapless 
approximation |2^. 
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FIG. 4. Radial eigen-functions R" (r) of the phonon-like 
basis states of the mean-field Hamiltonian operator 

versus radius. Depicted are a few representative 
states for the radial and angular quantum numbers ri^ — 1, 2, 10 
and 1 = 0,2. 



D. Time-dependent solution of the 
Hartree-Fock-Bogoliubov equations 

After studying some aspects of the stationary solutions 
of the generalized HFB equations Eqs. ( [zi]) and (72), 
such as the local densities, the eigen-states, or the energy 
spectra, we will now investigate the reversible real time 
evolution of the coupled condensate and noncondensate 
system, i. e.. 



^G> ^-iJ:G> +iG> 
at 



(75) 
(76) 



In contrast to the complete kinetic Eqs. ( p^ ) and (26), 
which account for coUisionally induced energy shifts and 
irreversible population transfer, Eqs. (|7|) and (|7^) con- 
tain only the first order reversible processes. As we will 
show in the following, this restriction implies constant 
occupation numbers P. 
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FIG. 5. Energy spectra versus radial and angular quan- 
tum numbers 1 < ti"^ < 8 and < Z < 6 for the phonon-like 
states of [l = 0, only), the particle-like states of Ea^ , as w/ell 
as the positive part of the energy spectrum of 11 (/ = 0, only) 
and HFB self-energy E. Dimensionless energies are measured in 
natural harmonic oscillator units. 



In Sec. IV A, we have shown that any admissible fluc- 
tuation matrix has to be of the form 



G> ^V+PVl + V_ (1 + P)V1, 



(77) 



where P represents the positive occupation numbers of 
the eigen-states V+ . This property is not only to be sat- 
isfied in equilibrium where the eigen-states coincide with 
the HFB states [see Eq.(74)], but in all instances. 

By formally integrating the reversible kinetic equa- 
tions, we can show that this structure of the fluctuation 
matrix is preserved at all times. Thus, Eqs. (^5|) and 
([76| ) define a consistent initial value problem. This sim- 
ple but important fact can be demonstrated easily by 
defining the formal solution in terms of a time-ordered 
exponential: 



T{t,to) =Texp[-i dt'T,{t')] 



(78) 



to 



It is obvious that the propagator matrix T{t,to) can be 
constructed in two steps by T{t,to) — T{t,ti)T{ti,to) 
since the self-energy is local in time (semi-group prop- 
erty). Moreover, it follows from the structure of the gen- 
erator that T{t, to) is a proper symplectic transformation 
CT3 — T{t,to) asT{t,toy at all times. Consequently, all 
occupation numbers of the general solution to the non- 
linear, initial value problem are constants of motion 



G>{t)^T{t,to)G>(to)T{t,to)^ 



(79) 



= V^t) P{to) Vlit) + V_{t) (1 + P{to)) Vlit), 
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with time-evolved basis states V{t) = T{t, to) V. Since 
T(t, to) represents a genuine symplectic transform, the 
eigen-states VHt) erg V{t) = a^V = (T3 remain or- 
thogonal. 

In the following figures, we illustrate these funda- 
mental facts that effectively define reversibility for any 
non-linear system in the generic form of Eqs. ( |75| ) and 
([TGI). In particular, we will use the self-consistent, fi- 
nite temperature solution for xi^ = 0) = ^-nd 
G^{t = 0) = G^{f3) as an initial value for the time- 
propagation at i = 0. It is obvious that this choice does 
not induce any change during the subsequent real-time 
evolution. However, ai t = 1, we suddenly distort this 
equilibrium solution by setting the anomalous compo- 
nent of G^{t = 1_) to zero, i.e., m{t = 1+) = and 
then propagate forward up to t = 4. It is important to 
note that the new G-^{t = 1+) is still a valid fluctuation 
matrix. 
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FIG. 6. Real time evolution of the total particle number 
= A'^'^ + A'' (solid), the number of particles in the con- 
densate A'^'^ (dash-dot), the noncondensate particles number A" 
(solid) and the trace over the anomalous fluctuations Tr{m} 
(dashed-dot). At t = 1, the equilibrium solution is suddenly dis- 
torted by setting fh{t = 1) — 0. Dimensionless time is measured 
in natural harmonic oscillator periods. 



In Fig. |6[ we have depicted the number of particles 
that occupy the condensate TV*^"^) = TT{f^'^\t)}, the non- 
condensate N — Tr{/(t)}, and the total particle number 
versus time. Time is measured in units 
of the harmonic oscillator period T = 2-k /uj. In contrast 
to these numbers that are genuine single-particle proper- 
ties, the anomalous fluctuations are a physical measure 
of the degree of two-particle correlations or squeezing 
p9|. For example, the total particle number fluctuations 
{{N — {N)Y) or, more specifically, the normally ordered 



density fluctuations (: /(x, x), /(y, y) :) would reflect the 
degree of squeezing of the quadrature components along 
certain directions x, y. While we have not evaluated such 
observables here, we have included the averaged strength 
of the anomalous fluctuations Trim} to represent their 
size. 

The most important feature in Fig. ^ is the exact con- 
servation of the total particle number during all phases 
of the evolution. The instantaneous change m. m{t = 1) 
does not affect it directly. But it can be seen that the 
relative partitioning of the particles between condensate 
and normal noncondensate is massively distorted by this 
sudden influx of energy and particles start to oscillate 
coherently between the coupled subsystems. 

In Fig. 0, we show the individual flrst order contribu- 
tions to the total system energy E = E'^'^^ -\- E[f]-\- E[rn\ as 
described in Eqs. (^6|)-(^8|). Again the most notable fea- 
ture is the exact conservation of total energy during the 
time-evolution. Due to the sudden change in the anoma- 
lous fluctuations a.i t — 1, the overall energy increases 
instantaneously by more than 400 Sw. 
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FIG. 7. Real time evolution of the total system energy 
E = E^"^ + E[f] + E[m] (solid), the energy of the particles 
in the condensate E'^'^^ (dashed-dot), the normal noncondensate 
energy E[f\ and the anomalous energy i?[m]. At t = 1, the equi- 
librium solution is suddenly distorted by setting m(t = 1) = 0. 
Dimensionless energy is measured in natural harmonic oscillator 
units. 



In Fig. ^ we show the complete spectrum of occupation 
numbers P(l < n'' < 14, < ^ < 6; i) that occur in 
the instantaneous fluctuation matrix G^(t), as defined 
by Eq. (^9|). From the logarithmic plot that covers 13 
decades, it can be seen that the occupation numbers are 
indeed numerically exact constants of motion. Due to the 
instantaneous change at i = 1, many of the high-energy 
occupation numbers that are quasi-degenerate before the 
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change, split up into a multitude of nondegenerate levels 
afterwards. Minute changes in the occupation numbers 
\6P\ < 10^^" at the end of the integration period identify 
the precision loss of the numerical calculation. 



10 



10 



10 



P 10° 



10 



10 



10 



P{n',l 



0.5 



1.5 



2.5 



3.5 



FIG. 8. Real time evolution of the instantaneous occupation 
numbers P(l < n'^ < 14, < Z < 6) that characterize the 
fluctuation matrix G^{t). Att = 1, (t) is suddenly distorted 
by setting fh{t = 1) = 0. 



The numerical results discussed in this section were 
obtained by discretizing the generalized GP Eq. (|7^) and 
the reversible part of the Boltzmann Eq. (^6|) with a 
standard finite element method |Q based on b-splines 
1^,^. The use of these finite-support, piecewise poly- 
nomial basis functions results in matrix representations 
of the kinetic and potential energy operators that are 
banded. Very efficient linear algebra algorithms can be 
employed in this case (LAPACK |5^]). In particular, we 
used a set of 400 b-splines on an equidistant radial grid 
from < r < 25 ~ 6 ttf • Eventually, we represented 
the condensate wave-function a (Fig. ^), the particle- 
like basis states |i) (Fig. ^), as well as the phonon-like 
basis sates |1^^^) (Fig. ^) in this b-sphne basis. Finally, a 
subset of quantum states was chosen (either |i) or 
with {(n'' > 1,1 > 0) : 2(n'' - 1) + I <= 18 ^ 2^} to 
evaluate the finite temperature sums or to propagate in 
real time. We have verified numerically that no particu- 
lar advantage can be obtained from either choice, as long 
as all of the relevant energies scales are resolved. The 
conservation of the total particle number (N) (Fig. |^), 
the total energy E (Fig. ^) as well as the instantaneous 
occupation numbers P (Fig. @) support this argument. 



V. THE IRREVERSIBLE EVOLUTION 



A. An ergodic equilibrium solution of the 
master-equation 

In the following discussion of the irreversible evolution 
of the kinetic equations, we will again try to elucidate the 
main physics by additional simplifications. In particu- 
lar, we will assume ergodicity for the normal fluctuations 



/i2 = ^12 /ei , and the anomalous fluctuations 



0. 



This physical approximation is appropriate for most ki- 
netic temperatures, except for a region close to T = 0. 
Within this limit, we are able to establish an important 
result for the stationary behavior of the condensed atomic 
gas, I.e., a canonical Bose-Einstein distribution for the 
noncondensate particles coexisting with an energetically 
lower-lying, coherent condensate mode. By imposing the 
restriction of vanishing anomalous fluctuations upon the 
kinetic equations, [Eq. (^ and (p^], we are left with 
the following equations of motion: 



—a = {-iU^ 



df' 



(80) 
(81) 



It is worth mentioning that in the real time evolution the 
rotating frame frequency fi is still an adjustable param- 
eter and not necessarily synonymous with the chemical 
potential. 

First, let us concentrate on the equation for the mean- 
field amplitude a . From Eq. ( |80| ) , we can see that the 
mean- field evolution consists of two distinct parts: a Her- 
mitian, number-conserving contribution ll^f, and a part 
that accounts for condensate number changing collisions 
out of and into the noncondensate, i.e., — T^. In 
stationarity, both processes have to vanish identically 



I a . 



(82) 
(83) 



//(!+/) (!+/)(! + /)/ • 

Given the constraint on the condensate particle num- 
ber, it is in principle straightforward to solve the Her- 
mitian eigen-value problem of Eq. (p2) and obtain an 
eigen- value fi. The later equation Eq. (|83|) poses a much 
more challenging constraint on the coupled system. In 
order to maintain a stationary state, the mean-field has 
to be "orthogonal" to all number changing processes. 
This means that the normal fluctuations have to adjust 
self-consistently with respect to the condensate wave- 
function. It is interesting to note that the solutions of 
this system are infinitely degenerate with respect to a 
global phase-rotation. In other words, the condensate's 
phase is not pinned down by any restoring force and is 
free to drift, consequently. However, the later, vector- 
valued condition of Eq. ( [8^ ) is satisfied identically if 

i4"l 2"l" j1"2"3"4" 



= 



(1 + /. 



,)(l + f\„)a,„ 

fe,,, 



' + £2" — — £4" 



(84) 
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vanishes for all components 1. Provided that energy con- 
servation is satisfied exactly, i. e., ey +£2" = /i-l-e4" , it is 
straightforward to verify that a canonical Bose-Einstein 
distribution 

- e.(e-l)-l ' 

with an inverse temperature /?, is the equilibrium solu- 
tion. In addition to this functional form of the distribu- 
tion that is dictated by detailed balance, it is required 
that all of the excitation energies are above the conden- 
sate energy, i. e., ei > ^. Thus the eigen-energy spectrum 
exhibits a finite gap. From Fig. |^ it can be seen that both 
the positive energy spectrum of the HFB operator E as 
well as the eigen-energies of E^r, are suitable candidates 
within this approximation. 



Second, the generalized Boltzmann equation, Eq. (81) 
is stationary if 







(^/7(i+/) + ^^/('•J /(!+/) 



+ r^-^-^,., )(! + /) (86) 



(i+/)(n-/)/ 



/<"'(! + /)/ (l + /)(l+/)/< 



Within the ergodic approximation, the Hermitian part 
is satisfied identically. On the other hand, there are 
two distinct types of coUisional relaxation processes in 
Eq. (p6[). There are the number-conserving in and out 
rates of the conventional quantum-Boltzmann equation, 
i-e., ^ff(i+f)(.'^ + f) - ^(i+f)(i+f)ff- Obviously, both 
rates match identically under the detailed balance condi- 
tions of Eq. (p5|) 



(l + /e,„)(l + /s,„)(l + /e,„)(l + /ej 



(87) 



Both of the two remaining distinct processes in Eq. (86) 



2r, 



(i + /)-2r 



f) f f, as well as the 



process rjj^(e) {1 + f) - T ^^^^j^f^^_^_j^ j(^) /, involve a con- 
densate particle in the pre-coUision or post-collision chan- 
nels. Thus, real particles will be transfered between 
the condensate and the noncondensate, until the rates 
are balanced. Analogous to the arguments that lead to 
Eq. (0), it can be shown that a canonical Bose-Einstein 
distribution is attained in equilibrium. 



VI. CONCLUSIONS AND OUTLOOK 

In this article, we have studied aspects of the reversible 
and irreversible evolution of a condensed atomic gas im- 
mersed in the noncondensate. By specializing the kinetic 
equations for a simple isotropic model, we were able 
to analyze the equilibrium solution, as well as the dy- 
namic nonequilibrium behavior numerically. In particu- 
lar, we obtained the excitation spectra of a finite tem- 
perature equilibrium. Moreover, we demonstrated the 



reversibility of the time-dependent HFB equations far 
from equilibrium and in the collisionless regime. This 
is tantamount to noting that the instantaneous occupa- 
tion numbers of the HFB modes are constants of motion. 
Finally, we studied the coUisional regime for an ergodic 
distribution of quasi-particles and showed that detailed 
balance is obtained in the full quantum kinetic theory 
with a self-consistent canonical Bose-Einstein distribu- 
tion. 

Based on this isotropic model, we can also obtain the 
collision rates that lead to a self-consistent equilibrium. 
However, such an analysis is still work in progress and 
results will be presented in future publications. 
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APPENDIX A: CANONICAL 
TRANSFORMATIONS 

A canonical transformation is an inhomogeneous linear 
combination of creation and destruction operators that 
preserves the commutation relation JlGf . In particular, if 
a and denotes a pair of Hermitian conjugated bosonic 
operators, such that 



\,dl] 



Si, 



(Al) 



then any affine linear transformation defines a new set of 
operators b and b by 



= T 



(A2) 



In an n-dimensional vector space, T represents a (2 n x 
2 n) dimensional matrix and d is a (2 n) dimensional vec- 
tor. Such a transformation is canonical if the new pair 
of operators also satisfies the commutation relation: 



[61,62] = Si,; 



(A3) 



More specifically, the transformation is unitary canonical 
if the new operators are H erm itian conju gate pairs, i. c., 
b ~ b^ . By inserting Eq. (A2) into Eq. (A3), one finds 
that the transformation matrices are a representation of 
the symplectic group Sp{2n): 



(A4) 



In addition, it can be shown that T* = aiT a\ and 
T^^ — (JsTVa. Here, we have introduced the (2n)- 
dimensional Pauli matrices a\ and as 
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1 

1 



(^3 



1 

-1 



(A5) 



APPENDIX B: CAUCHY-SCHWARTZ 
INEQUALITY 

For a positive semi-definite density operator a and a 
general operator L it follows that the expectation value 



{lV) = Tr{crLi:^} > 



(Bl) 



is never negative. Consequently, the co-variance matrix 
of Eq. (eJ) must be positive semi- definite u^G>u>^, 
as well. This can be easily seen, by considering a linear 
combination of two arbitrary operators A and i3, i.e., 
L = a A + f3B. By minimizing the positive expression 



Eq. (Bl), one obtains the Cauchy-Schwartz inequality as 
(ii^XBB^) > (Sit)(iB^). (B2) 



In particular, for the special choice of A 



and 



this implies that the magnitude of the 



B = al 

anomalous fluctuations is limited by 



(l + /ii)/22 > Imu? 



(B3) 
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